function I = SimpsonIntegration( fx, a, b, n )
f = inline(fx);
h = (b - a) / (2 * n);
I = f(a) + f(b);
x = a + h;
while (x < b)
    I = I + 4 * f(x);
    x = x + 2 * h;
end;
x = a + 2 * h;
while (x < b)
    I = I + 2 * f(x);
    x = x + 2 * h;
end;
I = h/3 * I;
end

